Reformulation of Smoothed Particle Hydrodynamics 
with Riemann Solver 1 

Shu-ichiro Inutsuka 

Department of Physics, Kyoto University, Kyoto, 606-8502, Japan 
E-mail: inutsuka@tap.scphys.kyoto-u.ac.jp 



Smoothed Particle Hydrodynamics is reformulated in terms of the con- 
volution of the original hydrodynamics equations, and the new evolution 
equations for the particles are derived. The same evolution equation of 
motion is also derived using a new action principle. The force acting on 
each particle is determined by solving the Riemann problem. The use of 
the Riemann Solver strengthens the method, making it accurate for the 
study of phenomena with strong shocks. The prescription for the variable 
smoothing length is shown. These techniques are implemented in strict 
conservation form. The results of a few test problems are also shown. 



1. INTRODUCTION 

Smoothed particle hydrodynamics (SPH, [14, 7]) is a fully Lagrangian particle method to describe 
the hydrodynamical phenomena. The Lagrangian particle methods are especially suited to hydrody- 
namical problems that have large empty regions and moving boundaries. Those problems naturally 
arise in engineering science as well as geophysics and astrophysics. A variety of astrophysical problems 
have been studied by SPH because of its simplicity in programming the two- and three-dimensional 
codes and its versatility of incorporating the various physical effects such as self-gravity, radiative cool- 
ing, and chemical reactions. A broad discussion of the method can be found in a review by Monaghan 
[18]. However, the inconsistency of the method is emphasized by Dilts [3, 4] who modified the method 
by means of the moving-least-square basis functions to obtain accurate derivatives regardless of the 
positions of the SPH particles. Another concern of the method is its poor description of the strong 
shocks. In the two- or three-dimensional calculation of colliding gases, particles often penetrate into 
the opposite side. This unphysical effect can be partially eliminated by the so-called XSPH prescrip- 
tion [17] that does not introduce the (required) additional dissipation but results in the additional 
dispersion of the waves. Therefore it is very desirable to construct a method that is simple and able 
to describe the strong shock phenomena accurately. 

In this paper, a new method for handling shocks in particle hydrodynamics is constructed. The 
force acting on each particle is determined by solving the Riemann problem (RP). This use of the so- 
called Riemann Solver is introduced as an simple analogy of the grid-based method [6]. The previous 
attempts to introduce Riemann Solver into the particle method failed to give the method an exact 
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conservation form [9, 10]. This paper describes how to include the exact Riemann Solver into the 
strictly conservative particle method. 

The alternative approach to including small but sufficient dissipation into the numerical solution is 
to find a good limiter to switch a dissipative method and a less-dissipative method [5, 21]. In principle, 
however, the switch is always an option to any numerical scheme including the present one, and its 
discussion is beyond the scope of this paper. 

Section 2 provides a description of the method where we derive the evolution equations for the 
particles in terms of the convolution of the original hydrodynamic equations with the so-called kernel 
function. The same evolution equations are derived from an action principle which is different from 
the previous ones [19, 20], which sheds light on the "hidden" approximation in the expression for 
the velocity field in SPH formalism. The detailed explanation for the implementation is described in 
Section 3. The usage of the Riemann Solver is analogous to the grid-based second-order Godunov 
scheme [24]. A variable smoothing length is also considered. Numerical examples involving strong 
shocks are presented in Section 4. Section 5 is for summary. 

2. THE METHOD 

We consider the following set of equations for non-radiating inviscid fluid: 

g (i) 
£ - >> 

du P „ 

P = ( 7 - l)pu, (4) 

where u denotes the specific internal energy, 7 denotes the ratio of specific heats, and other symbols 
have their usual meanings. 

In the standard SPH method, each particle has its own mass, and density at each particle's location 
is simply assigned by 

p i = ^2m j W(x i -x j ,h), (5) 

3 

where subscripts denote particle labels, rrij is the mass of the j-th particle, and W(x, h) is a spherically 
symmetric kernel function which is normalized to be unity if integrated in space: 

W(x,h) = W(-x,h), (6) 



/ 



W(x,h)dx = l, (7) 



and h is the parameter of the kernel function. For later convenience we also define the effective width 
h e ff of the kernel function W(x) by 

h 2 eS = 2 J x 2 W{x,h)dx. (8) 
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Although there are many possible forms for W(x, h), we use the following Gaussian kernel throughout 
in this paper: 



W(x,h) 



1 



-x 2 /h 2 



(9) 



In the above equation, d is the number of dimensions, and h is the so-called smoothing length. Note 
that h e g — h for this Gaussian kernel. We assume h is constant in space in Sections 2.1-3.2. The 
spatially variable smoothing length is considered in Section 3.3 and the subsequent Sections. 

The standard procedure to derive sets of evolution equations for SPH has been previously described 
[18]. In the following sections, we show how new evolution equations are derived from the original 
hydrodynamics equations. We then show that our new equation of motion can be derived from the 
action principle for an appropriately defined Lagrangian function. 

2.1. Key Concept 

In general any method of computational fluid dynamics inevitably brings errors into the solutions. 
Therefore the essential feature of the method can be characterized by how to introduce the errors into 
the solutions. For example, in many kinds of grid-based methods, the truncation of the Taylor-series 
expansion is the origin of the errors in the numerical solutions. In the SPH formalism we consider the 
convolution of the physical function f(x) with the kernel function, 



(f)(x) = J f(x')W(x' -x,h)dx'. 



(10) 



We use the symbol () to denote the convolution. Taylor-series expansion of / around x = Xi in Eq.(10) 
shows that the difference of (f)(xi) and f(xi) is the second-order in h c g: 



(f)( Xl ) = f( Xt ) + ^V\f( x ) + 0(hi s ). 



(11) 



Thus if we use (/) as an approximate solution for /, the errors of the second-order in h c g are introduced 
by the convolution with the symmetric function W. 
We also have 



<l^) = j dJ ^W { x'-x, h)d x> 

= - J f(x')-^W(x'-x,h)dx' 
J f(x')W(x' -x,h)dx', 



d_ 

dx 



(12) 



where we used the integration by part and Eq.(6). That is, the kernel convolution of V/ is just V(/). 
If we define the density distribution by the summation of the kernel functions at particle positions: 



p(x) = ^ rrijW(x — Xj, h), 



(13) 



then we have the following identities: 



(14) 



W(x - Xj, h) 
p{x) 



(15) 
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These equations are the key equations in the present method. Using Eq.(14) the value of the kernel 
convolution can be cast into the following expression: 



fi = (f)(Xi) = 



-I 



f(x')W(x' - x u h)dx' 
fix'), 

3 



mj W(x' - x, , h)W(x' - Xj , h)dx' 



(16) 



where we define 



in j 



/(*') 



p(x>) 



j-W{x' - x u h)W(x' - xj, h)dx'. 



(17) 



Note that mjA/jj is symmetric with respect to i and j. In this way the value (/j) of the physical 
variable / at the i-th particle position is expressed as the summation of the contributions (Afij) from 
the surrounding particles. This expression shows the spirit of the present method, although the actual 
detailed evolution equations are presented in the following sections. 
In contrast the standard SPH adopts the following form for A/jj : 



TO, 



p( Xj ) 



W(xi — xj,h). 



(18) 



This corresponds to the crude approximation, W(xj — x', h) ss 6(x' — Xj) in Eq.(17), where 5(x) is 
the Dirac's S funcion. The standard SPH also implicitly assumes the following approximate equation: 



(19) 



instead of Eq.(14), although this is very poor approximation. 

In the MLSPH scheme by Dilts, an equation analogous to Eq.(14) is provided by the moving- 
least-square basis functions, although he needed somewhat heavy operations to construct his basis 
functions. Note that our method does not restrict the functional form of the kernel function to satisfy 
Eqs.(13),(14). 

2.2. Equation of Motion 

To derive the appropriate evolution equation for particles we take the convolution of the equation 
of motion (EoM), Eq.(2), 



/ ^di^~ w<yX ~ x '' h ^ dx = ~ J ~pix) VP ^ w<yX ~ x ' j h>,dx ' 

The right-hand side of this equation can be transformed into the following expression: 

,p( x y 



(20) 




p(x) 
P(x) 



W(x - x', h) + \Vp(x)] W(x - x', h)\ dx 

p\x) J 



V W(x — x', h) — -Yp^- 



V m,jW(x — Xj , h) 



W(x-x',h) }dx 



= J2 m if {^W(x - x', h)W(x - Xj ,h) - W(x - x', h)VW(x - Xj ,h)} dx 



(21) 
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where we integrated by part and used the identity Eq.(14). 
Thus if we adopt the following equation for the evolution of the particle positions, 

/dv(x) f 1 

—^W(x - Xi ,h)dx = - J -^-jVP(x) W(x - x h h)dx, (22) 

we obtain the evolution equation that is consistent and spatially second-order accurate to the original 
hydrodynamical equation of motion: 



rriiXi = —rrii mj 



[ P(x) 


' d 


d ' 


J P 2 (x) 


dxi 


d Xj _ 



W(x-Xi,h)W(x-Xj,h)dx, (23) 



where the overdot means a time-derivative. The antisymmetric appearance of i and j on the right- 
hand side guarantees the conservation of the linear and angular momentum of the particle system. 
From this equation we can deduce some properties of the EoM of the particle system. For example, 
if the pressure distribution is constant in space, the acceleration vanishes exactly in Eq. (23), which 
is evident in Eq. (22). The evolution equation of MLSPH also has a similar property although the 
standard SPH does not [3]. 

In this way the main approximation introduced in the present method is simply expressed in Eq. 
(22), while the standard SPH formalism does not have such simple explanation. 

2.3. Energy Equation 

We can follow the evolution of barotropic fluid with P = P{p) only by the EoM derived in the 
previous section. However an additional evolution equation for energy is required to describe the 
evolution of non-barotropic fluid. The energy equation is also needed to handle shocks. The energy 
equation guarantees the conservation of the total energy of the system in the absence of other sources 
or losses of energy, such as radiative heating or cooling. On the other hand, the strict conservation 
property in the numerical calculation guarantees the accurate description of strong shocks. This is 
because the structure of shock wave is determined by the Rankinc-Hugoniot relation which is the 
direct consequence of the physical conservation laws. Thus we must derive the energy equation that 
has strict conservation property and a convenient form to include additional physical dissipation. 

As in Eq. (20), we multiply both sides of Eq. (3) by W(x — x', h), and integrate in space (with 
respect to aj), 

f ^l W (x - x', h)dx = - f ^4[V • v] W(x - x', h)dx. (24) 
J at J p(x) 

The right-hand side of this equation can be transformed into the following expression: 

- / ^t4[V • v] W(x - x 1 , h)dx 
J P(x) 

= - [ -4-r [V • Pv] W(x - x', h)dx + ( -!—[v ■ VP] W(x - x', h)dx 
J P\x) J p{x) 

(25) 

At this point, we again make an approximation according to the following equation (see also Eqs.(36),(38)): 
/ ~pjxj ^ ' VP ^ W<yX ~ Xi ' = J 'pJx) ][Xt ' W( - X ~ XU h ^ X + 

(26) 
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Then we have 



Ui = J W(x — Xi, h)dx 

« -y*-^y[V-Pu] W{x-x t ,h)dx + j -Lj[ Xi -VP] W{x-x t ,K)dx 

f \ W{x-Xj,h) ] 
= / P\v — Xi\-\/ dx 

J LP. 

= y^mj J [v — Xi] ■ [VW(x — Xi, h)W(x — Xj, h) 

—W(x - Xi, h)VW(x — Xj,h)] dx, 



(27) 



where we have again used the identities (14), (15). This form of the energy equation will be further 
modified to account for the physical dissipation in Section 3.2, and its conservation property will be 
shown in Section 3.4. 

2.4. Action Principle 

In this section we show another feature of the equation of motion in the present method. In the case 
without dissipation, the equation of motion for the particles in our method can be derived from an 
action principle. Therefore the time evolution of the particle system can be formulated in a canonical 
transformation, that might further enable the possible improvement of the time integration. 

For the moment, we consider barotropic fluid in which pressure is a function of density alone, 
P{p) = Kp 1 . In this case, the equation of motion (Eq.[2]) for the fluid can be derived by minimizing the 
action S which is the time integral of the Lagrangian function L expressed in Lagrangian coordinates: 



/ 



Ldt, (28) 



-fa--} 



pu 



dx. 



(29) 



If we use Eulerian coordinates, the appropriate Lagrangian density C must include constraint equations 
along with Lagrange multipliers [13, 22], which introduce additional complication. In contrast, the 
Lagrangian density shown above is simple and analogous to the classical particle mechanics form, owing 
to the use of Lagrangian coordinates. In addition we can derive the Hamiltonian and the evolution 
becomes the canonical transformation which is free from dissipation. Although this formalism is for 
the continuous media, we can use this formalism to derive the evolution equation for the system of 
particles. 

First we consider the density field expressed in the following equation: 



p(x) = ^ rrijW(x — Xj, h). 



(30) 



This continuous distribution of density is defined only by (the finite number of) the positions of 
particles, a;, (i = 1, 2, 3, N p ). Time-derivative of the above equation gives 



dp(x) 
dt 



d 



rrijXj ■ — — W(x — Xj, h) = — ^ m^Xj\IW{x - x ,h) = —V • F m , (31) 

A 3 A 
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where we defined 



F m (x) = ^ rrijXjW(x — Xj, h). 



(32) 



Comparing Eq. (31) with the ordinary continuity equation, we realize that F m defines the mass flux, 

F m (x) = p(x)v(x) (33) 
From this equation we can deduce the definition of velocity field given by the following equation: 

F m (x) _ Ej m jXj W{x - Xj ,h) 



Note that at x = x% the above equation gives 

v(xi) = Xi H ^ m,j[xj — Xi]W(xi — Xj,h). 



(34) 



(35) 



The difference of v(xi) and Xi is second order in h. In a one-dimensional case, for example, this can 
be shown by using an appropriate smooth function f x that satisfies f x (xi) = %i (i = 1, 2, N p ) as 



1 / df x 



Id 2 /* 



2 dx 2 



[xj - a; 4 ] 2 + 0([xj - .t,] 3 ) ^ W(a;i - xj, h) 



i fe 2 (df x d P , ia 2 /x 



Pi 2 [ 9x 9x 2 9a; 2 



p + 0(/i 4 ). 



(36) 



Eqs. (30) and (34) can be used for v(x) and p(x) in Eq. (29) to obtain the Lagrangian function 
that is defined only by the positions of particles. 




V - rrijXjW(x — Xj, h) 



2 V rrijW(x — Xj,h) 



K 

7^1 



^ rrijW(x — Xj , h) 



> dx. 



(37) 



This is the exact Lagrangian for the "fluid" of which density and velocity are constrained by Eqs. (30) 
and (34). 

Instead of adopting this complicated Lagrangian function, we now make an approximation for the 
velocity field by Taylor-series expansion of the velocity field: 

v(x) = v{x i ) + [x-x l ]-Vv{x i )^0{\x-x l \ 2 ) (38) 
= Xi + — y^ y mj[xj - Xi]W(xi - Xj,h) + [x - Xi] ■ Vv(xi) + 0(\x - Xi\ 2 ). 

From this we have 

J v(x)W{x - x u h)dx = x,+ 0(h 2 ). (39) 
With this equation we can simplify the kinetic term of Lagrangian function as follows: 

J ^pv 2 dx = ^ J F m - vdx 
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(40) 



Now we can define a simple Lagrangian function that is second-order accurate in h as 
L ncw = ^TO; ^xf - J u(x)W(x - Xi,h)dx . 



(41) 



We can derive the evolution equation for the particles by minimizing the action S = J L ncw dt . In 
fact the Euler-Lagrange's equation, 



d dL nev 
dt dx; 



dx. 



gives 



Xj 



f P(x) 


' d 


d ' 


1 P 2 (x) 


dxi 


dxj 



W(x — Xi, h)W(x — Xj 1 h)dx. 



(42) 



(43) 



The manipulation to derive Eq. (43) from Eq. (42) is explained in Appendix. This equation is the 
same as Eq. (23). 

This is the exact evolution equation for the system of which Lagrangian function is defined by 
Eq. (41). The space-symmetry of L ncw guarantees the conservation of linear momentum and angular 
momentum, which is also obvious in the anti-symmetric appearance of i and j in Eq. (43). 

If we start with more simplified Lagrangian function, 



in-,. 



£sph = 

i 

we would obtain the standard SPH equation [19], 



1.2 



(44) 



Xi 



3 



HI; 



„2 



dx, 



■W( Xi — Xj, h). 



(45) 



This Lagrangian function Lsph for the standard SPH can be obtained if we make a crude approxi- 
mation, W(x — Xi, h) « S(x — Xi) in Eq.(41). Thus the difference of the evolution equations, Eqs. 
(43) and (45) is related to the difference of the degrees of approximations to the original Lagrangian 
function of the fluid. In this context, the standard SPH is also assuming the approximation for the 
velocity field Eq.(39) as in the present method. 

Dilts [3] reported that the SPH approximation is derived by means of the Galcrkin approximation 
followed by the kernel approximation. In contrast to the Galerkin approximation, the above action 
principle has deep physical consequences such as the conservation of the linear and angular momen- 
tum of the particle system. In addition Hamiltonian structure of the method possibly gives further 
sophistication to the time-integration scheme (see, e.g, the symplectic integrator for the astronomical 
self-gravitating system [25]), although it is beyond the scope of this paper. 



3. IMPLEMENTATION 
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In this section we describe the numerical implementation of the method in detail. In Section 3.1 
we describe how to evaluate the integrals in the basic evolution equations. Section 3.2 is for the 
introduction of Riemann Solver into the method. Section 3.3 shows the prescription for the variable 
smoothing length. In Section 3.4 we show that the present method remains a fully conservative scheme 
even after the discretization in space and time. 

3.1. Convolution 

In this subsection we describe how to evaluate the spatial integrals in Eq.(23) and (27). First 
we have to know the distribution of p~ 2 (x) to calculate the integrals. That is, we have to construct 
appropriate interpolation (or extrapolation) function for p~ 2 (x) around each pair (i and j) of particles. 

For convenience, we define the s-axis which is along the vector Xi — Xj and has the origin at 
(xi — Xj)/2. We use s± to symbolically denote the components of the axes that are perpendicular 
to the s-axis, and define the s-coordinate system. The unit vector in the s-direction is ejj = (xi — 
Xj)/\xi — Xj\. We set Asij = S, — Sj = \xi — Xj\, where Si and Sj denote the s-components of the 
positions, Xi and Xj, respectively. 

We define the specific volume and its gradient as follows: 

v & = wr (46) 



VV{X) = -^) VP{X) = -^x) H m ^W(x Xj ). (47) 

We will make approximate function for V(x) in the s-coordinatcs. In the following subsections, two 
sets of equations with the different order of accuracy are shown. 

3.1.1. Linear Interpolation 

The most simple choice for the approximate function is the linear interpolation which is expressed 

as 

V{s) = ^=C^s + D i ^ (48) 



where 



_ V( Xi )-V( Xj ) 



_ vw + vw (49) 



From this we obtained 



v2 ( s ) = ^ = l C iJ s + D i,j} 2 > (5°) 



P 2 (s) 

With this interpolation, we can calculate the integral as follows: 



J p- 2 (x)W(x - Xi ,h)W{x - x 3 ,h)dx = V?j(h)W(xi - Xj ,V2h), (51) 



where we defined 



V 2 3 (h)^\h 2 Cl +Dl y (52) 



10 



SHU-ICHIRO INUTSUKA 



The above calculation corresponds to an approximation based on the linear expansion of 1/p 2 in the 
perpendicular direction to the vector x, t — Xj, that is, 

p- 2 (s) = V 2 (s) + s ± ■ W 2 (s) + 0(\s ± \ 2 ). (53) 

Note that the integral of [s± ■ VV 2 (s)] in Eq. (51) vanishes identically owing to the symmetry of the 
kernel. 

Next, we consider the following integral for arbitrary function f(x): 

[ ^\-W(x - Xi,h)W{x - Xj,h)dx. (54) 
J P \ x ) 

When we have some interpolation for the function f(x) and p~ 2 {x) , we can define the weighted 
average /,* • by the following equation: 

/ . W(x — x i ,h)W(x — x 1 ,h)dx^f*, / - W(x — Xi, h)W(x — Xj, h)dx. (55) 
J P A (x) ,J J p z (x) 

A straight-forward calculation of the above equation with the linear approximation for the function 
f(x) (w s[fi — fj]/Asij + [ft + fj]/2) and the specific choice (Eq.[50]) of the interpolation of p~ 2 (x) 
gives 

,* _ fi - fj * I fi+fj (r fi \ 

where 

2 y2 ■ l*>0 

That is, f*j is the value of the (linearly approximated) function / at the position s*j. 
The evolution equations now become the following: 

Xi = -2^m j P*V?j(h)-^W(x i -x j ,^h), (58) 



u t = -2^2 mj ([Pv]* -P*x i )V 2 j (h)-^W{x i -x j ,V2h), (59) 

3 1 

where the formal meanings of the P* and [Pv]* are the values at the position s*j of the linearly 
interpolated functions. In section 3.2, however, we will change the values of these quantities by 
considering the physical dissipation. 

3.1.2. Cubic Spline Interpolation 

Another convenient method for the approximation is based on the cubic spline interpolation of p~ x 
along the s-axis : 

V(s) = -iy = A itj s 3 + B^s 2 + C id s + D id , (60) 



where 



• ^ " (A S ,,)3 + (As . ,)2 . 
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2 As. 



C ^ 2 Asij l {Vi +V > h 



Dij = ^(V i + V j )--(V:-V')As i , j , 



1 . 1 

Vt = V( Xi ), 

Vj = V( Xj ), 

V( = eij-Wixi), 

V; = eij ■ VV( Xj ). (61) 



Then, 



y2{s) = 7(s) = ^ Ai ' js3 + B ^ + Ci ' jS + A ' j ] 2 ' (62) 

In this case Eq. (51) becomes 

p- 2 (x)W(x - Xi ,h)W(x - Xj, h)dx = V^(h)W(xi - Xj, V2h), (63) 

where we defined 



V ?M = l A h&A h + ^h\2A h3 a, + Bf d ) + + C&) + D\ y (64) 



64 hJ 16 

Eq. (57) becomes 



s i,j - 772 ' V 00 ) 



To avoid undershooting and overshooting of the interpolating function, we need to use linear interpo- 
lation (Eq.[48]) in the case V(V! < 0. 

3.2. The Usage of the Riemann Solver 

The description of shock waves requires (physical) dissipative process which is not considered in the 
fundamental equations (Eqs.[l]-[4]). The Godunov's scheme and its second-order sequel (MUSCL, [24]) 
and third-order sequel (PPM, [2]) use the exact Riemann Solver to include the (possibly) minimum 
and sufficient amount of dissipation into the method. At present these methods remain state-of-the- 
art grid-based methods for computational fluid dynamics. In this paper we introduce the exact and 
spatially second-order Riemann Solver into the particle method. 

The Godunov's scheme uses the result of Riemann problem (RP) at each cell interface in the 
calculation of numerical flux (see, e.g., [24]). Likewise, we want to use the result of RP at the vicinity 
of the middle point of the i-th particle and the j-th particle. This is achieved by modifying the values 
of P* and [Pv]* in Eqs. (58), (59). The finite difference expression for this is the following: 



A±i 

At 



= -I>^:,/^ 





' d 


d ' 


/ P 2 (x) 


dxi 


d Xj _ 



W(x — Xi, h)W(x — Xj, h)dx 



-2j2^P*V^(h)^W(x t -x h V2h), (66) 
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da;,; 



_9_ 

9^ 



W(x — Xj, /i)W(a; — ajj , h)dx 



d 



= 2j2rn j (P*v*-P*x*)V? J (h)—W(x i -x j ,V2h). 



(67) 



where A stands for the finite difference of each variable, and x* is the time-centered velocity of the 
i-th particle, and P*- and v* ■ are the results of RP between the i-th particle and the j-th particle. 
How to define and calculate these variables are explained below. 




FIG. 1. a) A schematic picture of the distribution of the function /(s) in the vicinity of the i-th particle and 
the j-th. The value of the weighted average /? . are determined by Eqs. (56), and its location s* . by Eq. (57) or (65). 
b) The setup of Riemann problem. We define the piecewise linear distribution of physical variable f(s) in both the i-th 
region and the j-th region. The initial values of each side of the one-dimensional Riemann problem are the average 
values in each domain of dependence. 

Figure (1) shows schematic picture of the distribution of the function /(s) in the vicinity of the i-th 
particle and the j-th. According to the original (grid-based) Godunov's scheme, it is natural to define 
the interface of the two regions for the RP at s*j. To make the spatially second-order method like 
MUSCL, we define the piecewise linear distribution of physical variable f(s) in both the i-th region 
and the j-th region. The gradient of f(s) in each region can be simply assigned by the gradient at 
each particle's position. The initial values on each side of one-dimensional Riemann problem is the 
average values of each domain of dependence, which is expressed as follows: 

s i,j ^s.i „ s i j 



PR = Pi + (I) 
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where 



v j = v j ■ ei,j- (69) 

and C Sy i denotes the sound speed at Xj. Thus we can solve the RP at the interface s*j. The detailed 
explanation for the solution of RP can be found in [24, 2], and we will not repeat it here. 

After solving the RP at s* •, we have the resulting pressure P*j and the projected velocity v*j. 
Three-dimensional (or two-dimensional) velocity t?* • is calculated as follows: 

V i,3 = V i,3 e h3 + i V *d ~ V i,3 e i,j}> ( 70 ) 

where 

Vij=Vi[l/2 + e]+Vj[l/2-e], (71) 



"t,3 



Vi[l/2 + e]+ Vj [l/2-e], (72) 
- - (73) 



These P*- and w* • are used in the calculation of the right-hand side of Eqs. (66) and (67). Note that 
Eqs.(71), (72), (73) are not needed in the actual calculation, because the bracketed second term on 
the right-hand side of Eq. (70) is perpendicular to eij and vanishes identically in Eq.(67). 

In the higher-order grid-based methods (e.g., MUSCL and PPM) we need to impose monotonicity 
constraint on the gradients of physical variable to obtain the stable description of the discontinuity. 
This is also true in the present method. Our experience shows that the monotonicity constraint only 
on the velocity field suffices to make the stable method. Thus in the actual calculation of our method, 
we mimic the monotonicity constraints by setting 

dv\ ( dv\ ( dv\ ( dv\ „ 

= if — x — < 0. (74) 



ds J i \ds J ■ \ds ) \ \ds / , 

The numerical scheme should be first-order in space at the surface of the shock wave. The shock 
surface tends to include a few particles in the actual calculation. Therefore, if the velocity difference 
of a certain particle pair corresponds to the sound speed divided by a small number, the pair is 
considered to be in the shock surface, and we need to use the first-order Riemann Solver for this pair. 
We implement this condition by setting 



ds J - \ ds 



if C s hock e it j ■ (vj - Vi) > mm(C Sti ,C s j), (75) 



3 



where / = p, P, v, and C s hock is a numerical constant corresponding to the number of particles at 
the shock surface. We adopt C s hock = 3 throughout in this paper. 

3.3. Variable Smoothing Length 

In the previous subsections, we assume that the smoothing length, h, is constant in space. In actual 
calculation we need to change h to enlarge the dynamic range of the spatial resolution. Therefore we 
have to extend the present method for the variable smoothing length. 

The formal derivation of the evolution equation with variable smoothing length is the same as 
Section 2.2, except that we should start with the definition of density as 

p(x) = TrijW(x — xj, h[x\). (76) 
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This definition of density corresponds to the so-called "gather" formulation of SPH [8] . The resulting 
evolution equations are the following: 

Sa = m j p *,j J ^) {^W(x - Xi , h[x])W(x - Xj ,h[x]) 

-W(x - Xi,h[x])VW(x - Xj,h[x])} dx, (77) 



E m j P iA v *J - v i) J {VW(* - Xi,h[x])W(x - Xj ,h[x]) 
—W(x — Xi, h[x])VW(x — Xj,h[x])} dx. 



(78) 



These equations are essentially the same as Eqs. (66) and (67), although the analytic integration 
for these equation is not possible even with the polynomial approximation for p~ 2 (x). Therefore we 
prefer to use simple approximation for the integral as in the following form: 



At 



= -y mj P* VfAh^—Wixi-Xj^hi) 

z — ' OXj 



+V* j (h j )—W(x i -x j ,V2h j 



(79) 



^ = -Y^rn j [P*v*-P*x;}\v^ j (h i )-^-W(x i -x j ,V2h i ) 

3 % 

+V* j (h j )-^ x -W(x i -x j ,V2h j ) 



(80) 



where, in spirit, we used hi for the half of the integration space which include Xj, and used hj for the 
other half. 

The above approximations assume that the h(x) should not vary much within the neighborhood 
of each particle. One possible way to determine the smoothing length with this constraint is the 
following: 



l/d 



(81) 



where 



Pi = E m j W ( X i - x JiK)i K = /ijCsmooth- 



(82) 



p* is more smooth than p itself if C smoo th > 1 • Numerical experiments shows that r) ~ 1 with 
Csmooth — 2 works fine in Section 4. The effective number of neighbors around each particle depends 
on the ratio of the smoothing length and the mean separation of particles at that particle. For 
example, we can usually ignore the contribution from the j-th particle to the i-th particle if their 
distance |x, — Xj\ is larger than 3hi, because exp(— 3 2 ) s=s 1.234 x 10~ 4 . Thus, in one-dimensional 
calculations, the number of neighbors for calculating Eq.(82) is 6?7C sm0 oth excluding the i-th particle 
itself. The number of neighbors becomes about 2877 2 Cg mooth in two-dimensional calculations, and 
about H3?7 3 Cg mooth in three dimensional calculations. 
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3.4. Conservation Property 

The final discretized form for the equation of motion for the i-th particle becomes the following: 



Ax, = -At^roj-f?. 



V&ihi^Wlxi-xj^hi) 



+V?A*i)-te:W{x i -x j ,y/2h j ). 

This calculation guarantees the conservation of the total momentum of the system 

TOiAi 4 = 0, 



(83) 



(84) 



because the terms in the square bracket in Eq. (83) are anti-symmetric in i and j. 
The final form for the equation of energy for the i-th particle becomes the following: 

A Ul = -A^m^*>* ;j -x*) V^ihJ-^-Wixi-xj^hi) 



+V* j (h j )—W(x i -x j ,V2h j ) 



where we defined the time-centered velocity of the i-th particle, 



(85) 



X^ X{ -\- ^A.x%. 

This expression guarantees the conservation of the total energy of the system, because 



1 .2 

2 X i + u '> 



Xi + Ax,] 2 + [m + Aui\ 



^ m l I A±i 



X{ ~\~ ^ ^X{ 



+ Au % 



+V? J (h j )-Q x -W(x i - Xj , V2hj) 



= 0. 



(86) 



(87) 



Thus the present scheme conserves energy exactly. This is in contrast to the ordinary energy equation 
of the standard SPH which is only accurate in the first-order in time (see, e.g., [1, 18]). We also note 
that the conservation of energy does not require constant smoothing length, which is obvious in Eq. 
(87). This is because the expression of the total energy in the numerical calculation ^ rrii(^x 2 + ui) 
does not explicitly depends on the choice of smoothing length. In other words, a sudden change of 
smoothing length at any timestep does not change the numerical value of the total energy. 

3.5. Overall Procedure 

In this section we summarize the actual procedures in the sequence of executions. The main loop 
of the time integration corresponds to Steps 1-4. 
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Step 0. Problem Setup 

We first setup the problem in the computer program. We appropriately place the particles to 
represent the density distribution that corresponds to the initial condition of the problem to solve. 
This may require some relaxation technique to find the appropriate positions of the particles [18]. We 
also determine pi, Vp,, and hi to start the main loop of the time integration. Various constants and 
the initial variables are calculated in this step. For example we determine the initial timestep At. 
Step 1. Gradient Calculation 

We calculate the gradient of physical variables P, v for the use in the Riemann Solver. 
Step 2. Source Term Calculation 

We calculate the RHS of Eqs.(83),(85) either by Eq.(52) or Eq.(64). The subroutine for the 
Riemann Solver is called once for every pair of particles to calculate P* and v* in Eqs.(83),(85). 
Step 3. Time Evolution 

We update Xi, x i} and Ui according to Eqs.(86),(83),(85) , 

Xi(t + At) = Xi(t)+x*At, 
±i(t + At) = Xi(t) + Axi, 

u t (t + At) = Ul (t) + A Ul . (88) 



Step 4. Density Updation 

According to the updated positions of particles, we update density distribution. The smoothing 
length of each particle is also updated. The timestep At for the next integration is also determined. 
We turn to Step 1 for the next time integration. 

4. NUMERICAL EXAMPLES 

The present method was tested on a variety of ID, 2D, and 3D problems, a few of which are 
described below. Other sets of test calculations will be described elsewhere. 

In determining the amount of the timestep At, we have to consider the Courant condition that is 
in spirit similar to the Courant condition for the grid-based Lagrangian methods, 



At = Ccfl min ■ 



17); 



Pi 




(89) 



The numerical experiments show that we can safely use Ccfl ~ 0.5 in most of hydrodynamical 
problems. Note that we don't need to consider the (effective) diffusion timescale that is related to the 
artificial viscosity adopted by the standard SPH and other particle methods. Thus At in the present 
method can be much larger than that in the other SPH methods. 

The following examples are calculated with the Fortran program in which the single precision real 
number is used. To accelerate the computation we use the data structure based on link lists [16]. 

Five cycles of iterations in the Riemann Solver is sufficient for the following test problems. 



4.1. Shock Tube 

In order to compare the capabilities of the present method and standard SPH we test our method 
against the shock tube problems, of which the exact solutions are available. The density variation in 
the first problem is not so large that we can use both the variable smoothing length and the constant 
smoothing length. The initial parameters of the problem are the following: 

PL = 1 , Pr = 0.5, 
P L = 1 , Pr = 0.2, 

v x ,l = , v XtK = 0, (90) 
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FIG. 2. Results of the shock tube problem in which the Mach number is 1.526. The right-hand side of the initial 
discontinuity includes 40 equal-mass particles. The open circles plot the snapshot at t = 0.2 by the present method with 
the cubic spline approximation for the convolution Eq.(64) and the variable smoothing length (rj = 1, C smoot h = 2.0). 
The solid lines correspond to the analytic solution. 



18 



SHU-ICHIRO INUTSUKA 



where the subscript "L" denotes the variable on the left-hand side of the initial discontinuity, and 
R denotes the right-hand side. The ratio of specific heats is 7 = 7/5. The Mach number of the 
resultant shock wave is 1.526. The value of the post-shock pressure is P* = 0.5099. Figure 2 plots the 
snapshot at t = 0.2 by the present method with the cubic spline approximation for the convolution 
Eq.(64) and the variable smoothing length (77 = 1, C smoo th = 2.0). The method with the linear 
approximation for the convolution Eq.(52) gives very similar result. Figure 3 plots the corresponding 
result of the standard SPH where we used the "standard" artificial viscosity (a = 1, /3 = 2) described 
in the review paper by Monaghan [18]. In both cases, the number of equal-mass particles is 80 (40) 
on the left(right)-hand side of initial discontinuity (i.e., Axl = 0.005, Axr — 0.01 ). The solid lines 
correspond to the analytic solution. In this problem, these two methods give similar results except 
for the contact discontinuity where the standard SPH produced a "wiggle" in pressure and specific 
internal energy. This is due to the inconsistency of EoM of the standard SPH (45). The present 
method produced the smooth distribution of internal energy at the contact discontinuity, and hence 
provide the almost constant pressure distribution at the density discontinuity. 
The total energy of the system (—0.4 < x < 0.4) is defined as 

,0.4 1 

/ -pu dx = 1.2. (91) 

The initial numerical value of the total energy of the particle system was the following: 

^rriiUi = 1.20000017 (92) 

i 

where the last few digits have no significant meaning in this Fortran single precision calculation. The 
final (t = 0.2) value of the total energy of the particle system was found to be 

+ Ui ) = 1.20003033 (93) 

i 

The relative error remains sufficiently small (AE/E < 10~ 4 ) when we change At in the other runs. 
Thus the error of the total energy is only due to the round-off error of the single precision calculation, 
and not due to the truncation errors in the numerical modeling of the evolution. This result guarantees 
the strict conservation property of the present scheme (see Section 3.4). 

The present method spent 0.17 seconds for the total 52 timesteps (At w 0.004) to compute with 
the Hewlett-Packard workstation C240 (PA-RISC 8200/236MHz), which corresponds to 3.27 x 10~ 3 
seconds per step and 2.72 x 10~ 5 seconds per particle. In contrast the standard SPH method used 956 
timesteps (At w 0.0002) for the stable evolution in this problem, because the artificial viscosity limits 
At. Indeed numerical experiments shows larger At causes unphysical oscillations in the solution. As 
a result the standard SPH method spent 1.35 seconds to finish the calculation. It corresponds to 
1.4 x 10~ 3 seconds per step and 1.2 x 10~ 5 seconds per particle. This means that the present method 
is about two times slower than the standard SPH for the required operations per particle but the 
present method tends to spend much smaller (total) computation time than the standard SPH does. 

Figure 4 plot the result of the present method where we used the constant smoothing length h = 0.01. 
The rarefaction wave on the left hand side is described less-accurately than in Fig. 2 as expected, 
because the smoothing length in high density region in Fig. 2 is smaller than the constant smoothing 
length in Fig. 4. The method with the constant smoothing length has no advantage over the method 
with the variable smoothing length even in this kind of shock tube problem where the density variation 
is not so large. 

Figure 5 plot the result of the present method with the first-order Riemann Solver where we set 
dp/ds = dP/ds = dv/ds = in the Riemann Solver. The sharp profiles at the shock front and the 
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head and tail of the rarefaction wave are smeared out as in the grid-based first-order Godunov method. 
In this calculation we don't need to calculate the gradients of P and v x . However the difference of the 
computation times with the first- and second-order method is negligible. There is no reason to use 
this first-order method. 




X 



FIG. 3. Results of the shock tube problem in which the Mach number is 1.526. The open circles plot the 

snapshot at t = 0.2 by the standard SPH. The solid lines correspond to the analytic solution. 
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FIG. 4. Same as Fig. 2 except that we used the constant smoothing length h = 0.01. 



SPH WITH RIEMANN SOLVER 



21 




-.4 -.2 .2 



X 



FIG. 5. Same as Fig. 2 except that we used the first-order Riemann Solver. 
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In Fig. 6, the results of the standard Sod's shock tube [23] are plotted. The initial parameters of 
the problem are the following: 

PL = 1 , PR = 0.125, 
Pl = 1 , Pr = 0.1, 

v x ,l = , v XiR = 0. (94) 

The ratio of specific heats is 7 = 7/5. The solid lines correspond to the analytic solution. The number 
of the equal-mass particles in the x-direction is 40 on the right-hand side of the initial discontinuity. 
We used the variable smoothing length (rj = 1, C smoo th = 2.0). 

4.2. Extreme Blast Wave 

To explore the capability of the present method, we test our method against an extremely strong 
shock tube problem. The initial parameters of the problem are the following: 

PL = 1 , PR = 1, 
P L = 3000 , P R = HT 7 , 
v x .l = , v x . R = 0, (95) 

where the initial pressure of the gas on the left-hand side is 3 x 10 10 times that of the right-hand side. 
The ratio of specific heats is 7 = 5/3. The Mach number is as large as 10 5 . We used 100 particles 
on each side of the initial discontinuity. Figure 7 plots the result of the present method where the 
cubic spline approximation for the convolution Eq.(64) and the variable smoothing length is used 
(77 = 1 , C sm0 oth = 2.0). Even in this severe problem, the present method gave stable and accurate 
results and is free from penetration problems. The pressure distribution shows a small wiggle at the 
contact surface. This is due to the approximation for the convolution. Figure 8 plots the result of 
the present method with the linear interpolation Eq.(52). The amplitude of the pressure wiggling at 
the contact surface is slightly larger than that of the cubic spline result. In order to eliminate this 
unphysical wiggling, we need to develop more accurate approximation for the numerical convolution 
than Eq.(64). Numerical experiments shows that the standard SPH cannot produce acceptable result 
at least with rj = 1, a w 1, and /3 w 2. The possible fine-tuning of the parameters in the standard 
SPH may enable the calculation in this case. Note that the present method can describe this extreme 
phenomena with the same parameters (77 = 1, C smoo th = 2.0) used in the other test problems, without 
any tuning of the method. 

The accurate calculation with the single precision real numbers in the Fortran program is due to the 
conservative formulation with the internal specific energy u instead of the total energy e (= v 2 /2 + u) 
or E (= pv 2 /2 + pu) which are usually used in the grid-based conservative numerical schemes (see 
Section 3.4). If we adopt e or E as the main variable for the integration, we have to do the subtraction 
to obtain u (= e — w 2 /2), which brings huge error into the numerical value of u in the extremely 
supersonic motion (i.e., when u <C e w v 2 /2). In this sense our method has potential advantage over 
the grid-based higher-order Godunov methods, at least, in describing extremely supersonic flows. 

4.3. Wind- Tunnel 

To test the present method against multi-dimensional problem, we adopt the following wind-tunnel 
problem [12], that can be in principle directly compared with the laboratory experiment. The geometry 
is a two-dimensional channel with a 15° wedge on the lower wall. A 15° expansion corner is also 
included. The inflow Mach number is two. This kind of problem is poorly suited for the particle 
method because a rigid-wall boundary condition must be set up at the wall of the tunnel. We present 
this test problem only to demonstrate the capability of our scheme. 
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FIG. 6. Results of the Sod's shock tube problem by the present method with the cubic spline approximation 

for the convolution Eq.(64) and the variable smoothing length (rj = 1, C smoot h = 2.0). The solid lines correspond to 
the analytic solution. 
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FIG. 7. Results of ID calculation for the extremely strong shock tube problem of which Mach number is about 
10 5 . We used 100 particles on each side of the initial discontinuity. The open circles plot the result of the present 
method where the cubic spline approximation for the convolution Eq.(64) and the variable smoothing length are used 
(rj = 1, C smoot h = 2.0). A small wiggle in the pressure distribution at the contact surface is due to the approximation 
for the convolution. 
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FIG. 8. Same as Fig. 7 except that we used the linear interpolation in the convolution Eq.(52). The amplitude 
of the pressure wiggling at the contact surface is slightly larger than that of the cubic spline result. 
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Fifteen Degree Wedge Channel Flow 




X 



FIG. 9. Results of calculation on "fifteen degree wedge channel flow" problem. Initially 32x96 particles are 

flowing inside the channel. Mach-number contours from 0.92 to 1.97 with an increment of 0.05 are plotted. 




FIG. 10. Same as Fig. 9 but positions of particles arc plotted. The crosses denote the positions of "ghost 

particles" that are introduced to mimic the rigid wall boundary condition. 
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FIG. 11. Results of calculation on "fifteen degree wedge channel flow" problem with the first-order Ricmann 

Solver. Initially 32x96 particles are flowing inside the channel. The contour levels are the same as Fig.9. 
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FIG. 12. Results of calculation on "fifteen degree wedge channel flow" problem with the standard SPH. The 

standard artificial viscosity with a = l,/3 = 2 arc used. Initially 32x96 particles are flowing inside the channel. The 
contour levels are the same as Fig.9. 
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FIG. 13. Results of calculation on "fifteen degree wedge channel flow" problem with the standard SPH. The 

standard artificial viscosity with a = 2, fi = 4 are used. Initially 32x96 particles arc flowing inside the channel. The 
contour levels are the same as Fig.9. 
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Figure 9 shows our method's result. Mach-number contours from 0.92 to 1.97 with an increment 
of 0.05 are plotted. Figure 10 shows particle positions. In this calculation, we realized the rigid-wall 
boundary condition by placing "ghost particles" in the wall. The inflowing particles outside of the 
left-hand side boundary were prepared at appropriate timesteps. The variable smoothing length is 
used (r/ = 1, C smoo th = 1-0). Initially, 32x96 particles were flowing inside the channel, which may be 
compared to the 32x96 grids finite difference calculation in [12]. Our results were satisfactory. For 
comparison, Figure 11 shows the result for the first-order Riemann Solver. The sharp features are 
smeared out. 

Figure 12 shows results of calculation with the standard SPH, where the standard artificial viscosity 
with a = 1, (3 = 2 are used. Note that the contour levels are the same as Fig. 9. Owing to the moderate 
strength of the shocks in this problem, the standard SPH could produce reasonable result except for 
the region just after the Mach stem on the upper wall. A small hollow semicircle seen at x ~ 0.5, 
y w 1 corresponds to the small value of Mach number (= 0.806) that is caused by an exceedingly 
large value of the post shock pressure. The calculations with larger values of the artificial viscosity 
parameters tend to remedy this pressure overshooting. Figure 13 shows the result of the standard 
SPH with a = 2, (3 = 4 where the minimum value of Mach number behind the Mach stem is 0.957. 
In this case, however, the sharp features are smeared out. Even in this two dimensional problem of 
moderate Mach number, the present method shows its potential ability in analyzing supersonic flows. 

It is clear that a grid-based method must be used in the actual study of this kind of steady rigid- wall 
boundary problem. This result is presented only to demonstrate the capability of this particle scheme. 
The present method becomes more useful when we study the hydrodynamical problems without rigid 
boundary that are common in astrophysics and space sciences. 



5. SUMMARY 

Smoothed Particle Hydrodynamics is reformulated by the formal convolution of the original hydro- 
dynamics equations, and by a new action principle, in which the second order (in h) approximation 
is used for the kinetic term of the Lagrangian function. The force acting on each particle is deter- 
mined by solving the Riemann problem for each particle pair. The prescription for variable smoothing 
length is also shown. These techniques are implemented in the strict conservation form. Numerical 
examples involving an extremely strong shock are shown. The other test calculations will be described 
elsewhere. 

Although the method with a spatially constant smoothing length is formulated in a rigorous manner, 
the method with the variable smoothing length is based on a crude approximation (Eqs. [79], [80]). 
A more refined technique for the variable smoothing length is to be studied. A better approximation 
than the cubic spline interpolation in the numerical convolution is required to eliminate the "wiggle" 
at the contact discontinuity in the blast wave problems (see Section 4.2). 

This paper presented a piece of concepts that are not discussed in detail in the literature. Those are 
the convolution of the original fluid equation (Eqs. [20], [24]), the definition and approximation for the 
velocity field (Eqs. [34], [39]), and the modification of the force due to dissipation (Eqs. [66], [67]). 
Incorporating these concept, the final evolution equation was cast into a form similar to the standard 
SPH equation. Therefore those concepts may enable the rigorous examination of the accuracy and 
stability of the method and may enable further modification. 

The numerical examples in Section 4 show that the present particle method based on Riemann 
Solver can handle severe problems with strong shocks, those might include the description of ex- 
plosion/implosion and supersonic jet phenomena. In this respect, further modification of the SPH 
method in modeling relativistic flows is promising with the help of relativistic Riemann Solver [15]. 
In adittion, the Lagrangian particle methods have advantage over Eulcrian grid-based methods, in 
describing chemically reacting (multi)fluid and radiatively heating/cooling fluid [11]. This is because 
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we can simply assign the chemical composition and entropy to each particle as a fluid element. This 
direction has a wide area of applications. 
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APPENDIX 

A.l. DERIVATION OF THE EQUATION OF MOTION 

In this Appendix we show the derivation of the equation of motion from the Euler-Lagrange's 
equation (Eq.[42]). 

= -^2 m k J §~ W ( X ~ Xk,h)dx 
f d 

—nii J u-^—W(x — x i} h)dx. (A.l) 
The first term on the RHS of this equation becomes the following: 

/du 
—W(x - x k , h)dx 

K 1 

/P dp 
f~faT W ( x - Xk > h ) dx 

K 

= -^m k f ^nii dW ( X — X%,h ^ W{x - x k ,h)dx. (A. 2) 

, J P OXi 



Before we manipulate the second term on the RHS of Eq. (A.l), we note that 

dW{x-x u h) _ dW(x-x u h) 
dxi dx 

We also obtain 



(A.3) 



,, .dW(x-Xi,h), [df{x). jr . ... .... 

dx — = ~ / dx ( X ~ Xi > h ^ dx > ^ A - 4 ) 

by the integration by parts and W(x) — > (\x\ — > oo) . 

Using the above identities, we can transform the second term on the RHS of Eq. (A.l) as 



f d 

~m,i / u— — W(x — Xi, h)dx 
J dxi 

f d 

+rrii / u— — W(x — Xi, h)dx 
J dx 

du, 



= —mi / ^-W(x — Xi, h)dx 
J dx 

= -mi J ^^W(x - Xi,h)dx 

= -m { [ rrij ~~~^~g~^ ^ W ( x ~ x i' h ) dx 

p . x 

= m * J ^ E m o dW %~ Xj,H) W(x - Xi , h)dx. (A.5) 
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As a result, Euler-Lagrange's equation gives the following: 



3 

3 J 



Xi, h)W(x — Xj,h)dx 



Xi,h)W(x — Xj, h)dx. 



(A.6) 
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